UNITS {
  (mV) = (millivolt)
  (mA) = (milliamp)
  (S) = (siemens)
}

NEURON {
  SUFFIX rgsoma
  USEION Na READ eNa WRITE iNa VALENCE 1
  USEION K READ eK WRITE iK VALENCE 1
  USEION NaP WRITE iNaP VALENCE 1
  NONSPECIFIC_CURRENT iLeak
  RANGE gKbar, gK
  RANGE gNabar, gNa
  RANGE gNaPbar, gNaP
  RANGE gLeakbar, eLeak
}

PARAMETER {
  gNabar = 0.030 (S/cm2)
  gNaPbar = 0.00025 (S/cm2)
  gKbar = 0.036 (S/cm2)
  gLeakbar = 0.0001 (S/cm2)
  eLeak = -64 (mV)
}

STATE {
  mNa  hNa  mNaP  hNaP  mK
}

ASSIGNED {
  v (mV)
  eNa (mV)
  iNa (mA/cm2)
  iNaP (mA/cm2)
  eK (mV)
  iK (mA/cm2)
  gNa (S/cm2)
  gNaP (S/cm2)
  gK (S/cm2)
  mNainf
  mNatau (ms)
  hNainf
  hNatau (ms)
  mNaPinf
  mNaPtau (ms)
  hNaPinf
  hNaPtau (ms)
  mKinf
  mKtau (ms)
  iLeak (mA/cm2)
}

BREAKPOINT {
  SOLVE states2 METHOD cnexp

  gNa = gNabar*mNa*mNa*mNa*hNa
  iNa = gNa*(v-eNa)
  gNaP = gNaPbar*mNaP*hNaP
  iNaP = gNaP*(v-eNa)
  gK = gKbar*mK*mK*mK*mK
  iK = gK*(v-eK)
  iLeak = gLeakbar*(v-eLeak)
}

INITIAL {
  rates3(v)

  mNa = mNainf
  hNa = hNainf
  mNaP = mNaPinf
  hNaP = hNaPinf
  mK = mKinf
  eNa = 55
  eK = -80
}

DERIVATIVE states2 {
  rates3(v)

  mNa' = (mNainf-mNa)/mNatau
  hNa' = (hNainf-hNa)/hNatau
  mNaP' = (mNaPinf-mNaP)/mNaPtau
  hNaP' = (hNaPinf-hNaP)/hNaPtau
  mK' = (mKinf-mK)/mKtau
}

UNITSOFF

PROCEDURE rates3(v(mV)) {
  TABLE mNainf, mNatau, hNainf, hNatau, mNaPinf, mNaPtau, hNaPinf,
   hNaPtau, mKtau, mKinf FROM -100 TO 100 WITH 32000

  mNainf = NaPvtrapB(-(v+35),7.8)
  mNatau = 1e-6
  hNainf = NaPvtrapB((v+55),7)
  hNatau = 30/(NaPeFun((v+50),15)+NaPeFun(-(v+50),16))
  mNaPinf = NaPvtrapB(-(v+47.1),3.1)
  mNaPtau = 1e-2
  hNaPinf = NaPvtrapB((v+59),8)
  hNaPtau = 2400/(NaPeFun((v+59),16)+NaPeFun(-(v+59),16))
  mKinf = NaPvtrapB(-(v+28),15)
  mKtau = 7/(NaPeFun((v+40),40)+NaPeFun(-(v+40),50))
}

FUNCTION NaPvtrapB(x,y) {
  LOCAL a
  a = x/y
  if (fabs(a)<1e-4) {
      NaPvtrapB = 1/(2+a)
  }else{
      NaPvtrapB = 1/(exp(a)+1)
  }
}

FUNCTION NaPeFun(x,y) {
  LOCAL a
  a= x/y
  if (fabs(a)<1e-4) {
      NaPeFun = 1+a
  }else{
      NaPeFun = exp(a)
  }
}

UNITSON